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ABSTRACT 

Magnetic fields play an important role in astrophysics on a wide variety of scales, rang- 
ing from the Sun and compact objects to galaxies and galaxy clusters. Here we discuss 
a novel implementation of ideal magnetohydrodynamics (MHD) in the moving mesh 
code AREPO which combines many of the advantages of Eulerian and Lagrangian 
methods in a single computational technique. The employed grid is defined as the 
Voronoi tessellation of a set of mesh-generating points which can move along with 
the flow, yielding an automatic adaptivity of the mesh and a substantial reduction 
of advection errors. Our scheme solves the MHD Riemann problem in the rest frame 
of the Voronoi interfaces using the HLLD Riemann solver. To satisfy the divergence 
constraint of the magnetic field in multiple dimensions, the Dedner divergence clean- 
ing method is applied. In a set of standard test problems we show that the new code 
produces accurate results, and that the divergence of the magnetic field is kept suffi- 
ciently small to closely preserve the correct physical solution. We also apply the code 
to two first application problems, namely supersonic MHD turbulence and the spher- 
ical collapse of a magnetized cloud. We verify that the code is able to handle both 
problems well, demonstrating the applicability of this MHD version of AREPO to a 
wide range of problems in astrophysics. 

Key words: methods: numerical, magnetohydrodynamics, turbulence, 
stars: formation 



1 INTRODUCTION 

In many astrophysical systems, the gas is partially or fully 
ionized, such that its conductivity is very high while its in- 
ternal viscosity is very low, implying that it can be treated in 
the limit of ideal magnetohydrodynamics. This is for exam- 
ple the case for the diffuse gas found in clusters of galaxies or 
in the halos of ordinary galaxies, and in most of the volume 
of the interstellar medium, where magnetic fields probably 
play an important role in regulating star formation. Also, 
magnetic fields are thought to be a key ingredient in mediat- 
ing angular momentum in accretion disks through the mag- 



netorotational instability ( [Balbus Hawley 1998| . Magnetic 
fields are also critical for launching relativistic outflows in 
some compact objects, and likely have a major influence on 
core-collapse supernovae. 

There is hence ample motivation to outfit numerical 
codes for hydrodynamics with an additional treatment of 
magnetic fields, prompting the development of a large num- 
ber of such codes in astrophysics ove r the years (e.g. [Stone k, 



Norman|1992l|Ziegler Yorke| 1997) [Brandenburg Dobler 



2002'; Trom ang et al.|[2006| jRosswog Price[ [20071 [Stone 
et al. 2008 : jPolag Stasyszyn|[2009| [Collins et al.|[2010[ ). 

Unfortunately, an important technical complication makes 



this far from straightforward. While the continuum equa- 
tions of ideal MHD preserve V • B = in an initially diver- 
gence free field, this is not necessarily the case for discretized 
versions of the equations. Here the relevant difference oper- 
ators will tend to pick up or produce a non-vanishing value 
for V • B, and once such a spurious 'magnetic monopole' 
has been produced, it has the tendency to become quickly 
larger in any non-trivial MHD flow, rapidly rendering the 
calculated solution unphysical. Simply ignoring this prob- 
lem is not a viable strategy, except for a small class of very 
simple and well behaved test problems. 

Much work has therefore been done on developing dis- 
cretization schemes for MHD that circumvent this problem. 
A detailed comparison of different approaches can be found 
in [T6th[ ([2000 ) . All approaches either try to construct the 
discretization such that it inherently forces V • B to be zero 
or modify the equations of ideal MHD to keep it small. 

A straight-forward way to get rid of V • B field compo- 
nents has been proposed by Brackbill & Barnes ( 19801). By 
projecting out the V • B component through a Helmholtz 
decomposition or a similar field cleaning technique one can 
regain a divergence-free magnetic field at any time. 

Another approach is to add additional terms that diffuse 
away the error once it appears (see, e.g. Powell et al.[[1999 
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Dedner et al.||2002 Keppens et al.|[2003 ). In particular the 
so-called Dedner cleaning technique ([Dedner et al.|2002| has 
proven to be quite robust and in most cases yields results 
that are of comparable accuracy as those obtained with the 



"constrained transport" scheme of |Evans & Hawley (1988) 
which fulfills the divergence constraint by construction. Here 
the transport of magnetic flux at the discretized level of the 
equations is carefully formulated in terms of loop integrals 
over electric fields at the mesh-edges, such that the magnetic 
flux and the V • B = constraint are conserved to machine 
precision. 

Unfortunately, the constrained transport method is only 
easily tractable for Cartesian meshes, and it is presently un- 
clear whether it can be adapted to dynamic unstructured 
meshes such as the one we consider in this paper. We shall 
hence resort to the Dedner cleaning technique. We note that 
yet another possibility for working around the V • B prob- 
lems lies in using the vector potential (e.g. Price||2010 ), or 
the so-called Euler potentials (e.g. Rosswog & Price 2007| 
Dolag &: Stasyszyii|[2Q09 ) . While both of these approaches 
formally guarantee V • B = 0, they are met with a number 
of serious problems in practical applications, in particular 
with respect to correctly accounting for magnetic dissipa- 



tion in turbulent flows (Brandenburg 2010), hence we do 
not consider them here. 



The AREPO code introduced by |Springel| ( |2010a| ) rep- 
resents a novel type of astrophysical simulation code, taking 
on an intermediate role between traditional Eulerian mesh- 
based hydrodynamics and mesh-free smoothed particle hy- 
drodynamics. In AREPO, a Voronoi tessellation is employed 
to construct an unstructured mesh, which is then used to 
solve the equations of ideal hydrodynamics with a second- 
order accurate finite volume scheme based on Godunov's 
method. A Cartesian mesh is in fact a special case of a 
Voronoi tessellation, and here AREPO 's basic fluid dynami- 
cal approach is equivalent to the widely employed MUSCL- 
Hancock scheme ( van Leer||1984| Toro||1997| ). However, one 
crucial difference is that in AREPO the mesh-generating 
points can be moved arbitrarily. In particular, they can be 
moved with the local flow velocity such that the mesh fol- 
lows the motion of the gas, providing quasi-Lagrangian be- 
haviour. In this mode, advection errors are greatly reduced 
compared to ordinary Eulerian techniques while at the same 
time their accuracy for shock waves and fluid instabilities is 
retained. This approach is then also not limited by global 
timestep contraints coming from the largest velocity in a 
system; instead, only the local timestep constraints apply, 
making the scheme ideal for systems with large bulk veloci- 
ties, such as accretion disks, or galaxy collisions. 

We note that these attractive features have already led 
other groups to start developing new codes based on simi- 
lar design principles. For example, TESS by Duff ell Mac-| 
[F^en| ( |20lT] ) is a new moving-mesh code that includes an 
extension to relativistic hydrodynamics and magnetic fields, 
albeit so far only in 2D, in serial, and without a control of 
the V • B error. 

In this study, we present the details of our new imple- 
mentation of MHD in the massively parallel AREPO code. 
We focus on a concise description of the numerical imple- 
mentation of the MHD part of the code, and a discussion 
of a set of basic test problems. A detailed account of the 



mesh construction algorithms and parallelization strategies 
can be found in Springel (2010a). 

The paper is structured as follows. Section [2] describes 
our implementation of MHD in the AREPO code. It is fol- 
lowed by a discussion of several standard test problems for 
MHD in Section |3] In Section [4] we show that the code can 
be used to simulate MHD turbulence. Additionally we apply 
the code to an important prototypical astrophysical applica- 
tion, the collapse of a magnetized cloud, in Sectionjs] Finally 
we give our conclusions in Section |6] 



2 IMPLEMENTATION 

2.1 The equations of magnetohydrodynamics 

The equations of ideal magneto-hydrodynamics can be writ- 
ten as a system of conservation laws, 

dU 

dt 

for a vector of conserved variables U and a flux function 
F(U), which are given in the local restframe by 



+ V-F = 0, 



(1) 



U : 



/ P \ 

pv 

pe 



F(U) 



pvv^ + p 
pew + pv — 



-BB^ 
B (v B) 



Bv' -vB' 



(2) 

Here, p — pgas + |B^ is the total gas pressure and e = 
u -\- + ^B^ is the total energy per unit mass, with u 
denoting the thermal energy per unit mass, p, v and B give 
the local gas density, velocity and magnetic field strength, 
respectively. 

For B = 0, these equations reduce to ideal hydrody- 
namics, which is treated by our version of AREPO in the 
same manner as described by Springe ij ([2010a ) , in particular 
with respect to the technical details of Voronoi mesh con- 
struction, gradient estimation, and parallelization. In the 
interest of brevity, we will therefore restrict ourselves in the 
following to a description of the aspects relevant for the new 
MHD implementation. 



2.2 Solving the Riemann problem 

To estimate the fluxes over an interface we follow Godunov's 
approach and solve an approximate Riemann problem nor- 
mal to the interface. We obtain the initial left and right 
state at the interface by spatially extrapolating the prim- 
itive variables at the centres of both neighboring cells to 
the mid-point of the interface, and by predicting them half 
a timestep forward in time. While the time extrapolation 
is done as described in |Springel| ( |2010a| ) , for the extrapola- 
tion in space we use a refined approach based on |Darwish fc| 
Moukalled ( 2003 ) . To extrapolate a primitive variable from 



the left cell to the interface, according to their proposal one 
first defines the scalar r by 



2V(/)L • TLR , 

= — : 1, 



(3) 



where (/)l,r are the values of the primitive variable (j) at the 
centers of the cells left and right from the interface, flr is 
the vector from the center of the left cell to the center of the 
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right cell and V0l the gradient of (j) at the center of the left 
cell. The value 0l of (j) on the interface extrapolated from 
the left cell can then be calculated as 



/^L + -^(rL) 



0. 



(4) 



where ^ (rL) is the slope limiting function. The extrapola- 
tion from the right cell to the interface is done likewise. Our 
usual choice for ^ is the van-Leer limiter ( van Lee"r||1974 ), 
defined by 



r + \r\ 
l + |r| 



(5) 



but other choices are in principle possible as well. Note that 
this slope limiting procedure recovers the common formula- 
tion for a Cartesian grid as a special case. 

To solve the Riemann problem in the MHD case, we use 
a three-step approach that is based on the use of up to three 
different Riemann solvers in the following order: 



1) The HLLD solve r ([Miyoshi k Kusario|[2005 ) 

2) The HLL solver ( [Harten et al.|1983l ) 

3) The Rusanov solver ( |Rusanov|l961| ) 

Each of these Riemann solvers returns fluxes over the 
interface as well as the state of the primitive variables at the 
interface. Starting with the HLLD solver, we check whether 
the pressure p* of the solution at the interface is positive. If 
this is not the case, the solution is not valid and is discarded. 
Instead, we revert to the HLL solver and try again. If p* is 
still negative, we also discard this solution and finally employ 
the Rusanov solver which guarantees a valid solution. We 
found this approach to provide a good compromise in terms 
of accuracy and robustness. 



2.3 The divergence constraint 

In addition to satisfying the MHD equations that govern the 
time evolution of the magnetic field, the field always has to 
obey the constraint 



V - B = 0. 



(6) 



It can be shown that the differential form of the MHD evo- 
lution equations keep a magnetic field divergence free if it 
has been in this state initially. Numerically, however, due to 
discretization errors this is not the case in general. In one 
dimension, the constraint reduces to dBx/dx — 0, which 
is equivalent to requiring a constant magnetic field in the 
x-direction, and is hence easily fulfilled. In two- or three di- 
mensions, however, it is in general non-trivial to satisfy the 
divergence constraint. 

Today, the most widely used approach to ensure V • B = 
in numerical MHD is the constrained transport scheme 
( Evans Hawley 1988 Gardiner Stone|20 05 ) which guar- 
antees that the divergence of the magnetic field remains zero 
by construction. However, it is very difficult, if at all possi- 
ble, to adapt this scheme to an unstructured moving grid. 
Therefore we here use an alternative approach proposed by 
Dedner et al. (2002), with the goal to keep the divergence of 



the magnetic field always sufficiently small rather than guar- 
anteeing it to be exactly zero. This can be viewed as the next 
best alternative compared to constrained transport. 



We apply the GLM-MHD approach ( [Dedner et al.|2002 ) 



by introducing an additional (conserved) scalar -0 which is 
related to the divergence of the magnetic field. The new vec- 
tor of conserved quantities and the new fiux function include 
-0, and are given by 



U : 



/ P \ 
pv 

pe 

B 

V ^ / 



F(U) 



pw'^ + j 
pev + pv - 



pv 



\ 



-BB^ 
B (v B) 

Bv^ - vB^ + tpl 



\ 



J 



(8) 



(7) 

Here Ch is a positive constant. In addition, we adopt a source 
term leading to an exponential decay of ip, 

d±^_cl 
dt 

where a second constant Cp has been added. We include this 
source term in the time integration using operator splitting. 
As shown in Dedner et al.| (|2002), it is possible to solve the 
equations for the fiuxes of Bx and separately, and to use 
an ordinary MHD Riemann solver in a second step applying 
the value of Bx at the interface that is obtained in the first 
step. The fiuxes are then formally given by 



Bx 



(9) 



There are no physical constraints on the choice of the 
two constants Ch and Cp but they have to be chosen such 
that the scheme is stable and prevents the divergence of the 
magnetic field from becoming so large that the solution of 
the problem is affected. For most standard test cases, both 
goals can be fulfilled without much of a problem. For turbu- 
lent motions, however, this is more challenging. In particu- 
lar, the choice of c^, which acts as the velocity with which 
the divergence of the magnetic field is advected away, re- 
quires some care. If it is too small, the transport will not be 
efficient enough and the scheme can become unstable due to 
a local build up of a large divergence error in the magnetic 
field. Conversely, if it is set too large, the physical solution 
of the problem can be impacted. 

A further difficulty is that the AREPO code evolves all 
cells on individual timesteps according to their local time- 



stepping criterion, as described by Springel (2010a). Equa- 
tion (|9| directly shows that Ch must be chosen globally, be- 
cause a different Ch at two opposite interfaces of a cell would 
lead to a net fiux of ^ even for a constant homogenous mag- 
netic field. The same argument demands that we may change 
Ch only when all cells are synchronized. As we do not know 
at that point in time how large it needs to be to keep the 
scheme stable until the next update opportunity, we assign 
the maximum signal velocity of all cells to Ch- Assuming that 
our grid motion is close to Lagrangian, the maximum signal 
velocity is given by the fastest magneto-acoustic wave 



Ch — maxi (c/) . 



(10) 



Since Ch acts as a velocity and will be larger than the local 
signal velocity for many cells, we need to introduce an ad- 
ditional limitation to the timestep of a cell, similar to the 
Courant-criterion, to keep the ordinary algorithm stable: 



AU < Cc 



Ch 



(11) 



Here, Ccfl is the same constant used for limiting the hy- 
drodynamical timestep and is the effective radius of the 
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cell. The choice of Cp turns out to be less problematic. As 
proposed by Dedner et al. (2002), we simply adopt 



Cp = V0.18 Ch. (12) 



2.4 MHD on a moving grid 

The fluxes described above are valid only for a static grid. 
If the grid itself moves, additional terms will be needed to 
account for the movement. The flux over an interface moving 
with velocity w can be described as a combination of the flux 
over a static interface and an advection step owing to the 
movement of the interface: 



F™(U) 



F.(U) 



pv 

pvv^ + p 
pey + pv - 



V 



-BB' 
B (v B) 

Bv^ - vB^ + V^/ 



J 



pvw 
pew 
Bw^ 

'0W 



(13) 

However, the straightforward approach to use the approx- 
imate solution of the Riemann problem in the rest-frame, 
Fs, and then to advect the state on the interface with the 
velocity of the interface, turns out to be unstable. This is 
mainly a result of our usual choice of the grid velocities, 
which cause the velocity of the interface to be very close to 
the fluid velocity at the interface, w v, and the mass flux 
over the moving interface is accordingly close to zero. As the 
mass flux over the interface in the rest frame is calculated 
from an approximative Riemann solver, there is always a 
small error. And because the mass flux in the moving frame 
is very small, this error can cause a sign change, destroying 
the upwind property of the scheme and making it unsta- 
ble. This behavior can be avoided by calculating the fluxes 
in the rest-frame of the interface. There the new vectors of 
conserved variables and fluxes are 



u' 



p (v — w) 
pe' 
B 



\ 



(14) 



F'(UO 



p(v- 
p (v — w) (v — w) 
pe' (v — w) + p (v — w' 
B(v-v 



-B((v- w) 
w) B^ + ^pi 



B) 



with e 



\ clB 

Q2 
Q3 
Q4 

V / 

(15) 

— + I (v — w) .To get the fluxes in the rest 



have to add some additional terms: 
F„(U) = F,(U) - Uw"^ = F'(U') 

/ 

pw (v — w) 



+ 



V 



p (vw) (v — w) — ^w^ (v — w) + pw — B (w • B) 
-wB^ 




/ 

(16) 

Finally, we can resubstitute the fluxes obtained from solving 
the Riemann problem in the moving frame into the addi- 
tional terms: 



F„(U) = F3(U)-Uw'' 
/ 

= F'(U') + 





wQf 

-wB'' 



V 







(17) 



With this procedure, the discretized MHD equations always 
retain their upwind character, and a stable evolution on a 
dynamic Voronoi mesh is obtained. 



3 TEST PROBLEMS 

In this section we study a number of standard test problems 
that are essential for validating the accuracy of a numer- 
ical MHD implementation. We consider first simple mag- 
netic shock tubes and then move to more demanding two- 
dimensional problems, each designed to test different aspects 
of the code. 



3.1 Shocktube 

The results of a standard one-dimensional shock-tube prob- 



lem with the initial conditions described by Keppens ( 2004 ) 



are shown in Fig.[T] This shock-tube problem has the advan- 
tage that all seven MHD waves are present. For comparison 
it is run for static and moving grids with two different Rie- 
mann solvers each (HLLD and Lax-Friedrich). Overall, all 
runs show good agreement. Upon close inspection, however, 
there are some differences. As expected, the HLLD-solver is 
less diffusive than the Lax-Friedrich solver. In addition, the 
moving grid behaves particularily better at the slow shock 
and the contact discontinuity, similar to previous studies 
van Dam Zegeling||2006 ) 



frame of the mesh from the fluxes in the moving frame we 



We note that in contrast to this one-dimensional shock- 
tube, all multi-dimensional problems we discuss next are run 
with an activated divergence cleaning scheme. Also, they all 
apply periodic boundary conditions and use 7 = 5/3 for the 
gas. 



3.2 Advection of a magnetic loop 

In this test problem, a magnetic fleld loop is advected by a 
constant velocity fleld. The magnetic fleld itself is too small 
to be dynamically important. The initial conditions are de- 
flned in a box of extension 0<x<2,0<2/<l and are 
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Figure 1. Results of a one-dimensional shock-tube described by |Keppens| ( |2004| with 7 = 5/3. The left and right initial states are given 
by {p^p^Vx^Vy^Vz^Bx^By^Bz) = (0.5,1.0,0,1.0,0.1,1.0,2.5,0) and (0.1,0,0,0,0,1.0,2.0,0), respectively. Shown are numerical solutions 
at t = 0.08 for a resolution of 250 cells for p (top row, top right panel zooms into the density discontinuity), Vx (middle left), Vz (middle 
right). By (bottom left), and p/ p (bottom right). The plots show results for runs with a static/moving grid using the HLLD Riemann 
solver (green dotted/black dashed line) and the Lax-Friedrich solver (red/blue dots). The position of the dots is given by the centers of 
the cells. 



given by 



V -■ 

V : 
A : 



(sin(7r/3), cos(7r/3), 0) 
(0,0,max(0.001 x (0.3 -r),0)) 



(18) 



where r is the radial distance to the center of the loop and 
A the vector potential. The initial magnetic field is directly 
calculated from the vector potential. We use a hexagonal 
grid with 1600 x 800 cells in total. 

In Figure [2] we show the magnetic field strength and its 
divergence error after the loop has crossed the box several 
times. The loop is preserved extremely well, and advection 
errors are very small, which highlights the principal advan- 
tage of the moving mesh code. Notice that the logarithmic 
plot of the magnetic energy density in Fig. [2] demonstrates 
that the smearing at the edge of the loop is quite small. 
More importantly, the relative error of the divergence of the 
magnetic field is smaller than ^ 10~^ for all parts of the 
magnetic loop except the very outer edge where the mag- 
netic field drops to zero. The error at the center of the loop 
is caused by insufficient angular resolution there. 



For a quantitative check of the advection properties of 
our MHD code we compare in Fig.[3]the time evolution of the 
magnetic energy with simulations on a static grid with the 
same resolution and a higher resolution by a factor of two. 
The moving mesh implementation conserves the magnetic 
energy significantly better than a simulation on a static grid 
of the same resolution. In fact, it is even a little bit better 
than a simulation on a static grid with a twice better spatial 
resolution. 



3.3 Magnetic blast wave 

The magnetic blast wave problem consists of an initially 
circular, overpressurized region in a magnetized fluid. The 
initial conditions are given by 



10 r < 0.1 
0.1 r < 0.1 



(0,0,0) 



(19) 
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Figure 2. Magnetic energy evolution of a field loop advected 
over a moving grid with a resolution of 1600 x 800 cells. The two 
panels on top show the magnetic energy density at times t = 
(first row) and t = 18 (second row) on a linear color scale, while 
the third panel repeats the t = 18 result with a logarithmic color 
scale. Finally, the bottom panel shows the relative importance of 
the divergence error of the magnetic field by comparing it to the 
magnetic field itself for all cells with |B| > 10~^. Here, d is the 
approximative size of a cell calculated from its volume assuming 
that it is spherical. 
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Figure 3. Time evolution of the magnetic energy in an advected 
field loop. The change of the magnetic energy is shown for a 
moving grid of 800 x 400 cells (black), a static grid of 800 x 400 
cells (red) and a static grid of 1600 x 800 cells (blue). 
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Figure 4. Linear density maps of a two-dimensional spherical 
blast wave at times t = 0.2 (left) and t = 3.75 (right). The density 
mapped to the color scale ranges from 0.1 to 3.0 (left) and 0.1 to 
2.4 (right). The system is resolved with 512 x 768 cells. 



3.4 Orszag-Tang vortex 



in a region of size 0<x<l, 0<2/< 1-5. Fig. [4] shows 
the density soon after the start of a simulation with a res- 
olution of 512 x 768 cells, as well at a later time when the 
shock has already interacted with itself after crossing the 
box. The early snapshot shows how the shock wave becomes 
anisotropic as the magnetic field introduces a preferred di- 
rection. 

The code captures the complex shock dynamics very 
well and manages to resolve quite sharp interfaces, such that 
the results favourably compare to corresponding calculations 
with other MHD cod es (see, e.g.|Fromang et al.|2006| |Stone| 
et al.|2008|[Rosswog Price|2007t . Note that there are some 
corrugations in the red regions behind the shock front in left 
panel of Fig. ^ that originate from small distortions in the 
grid. This can be cured by adding a refinement scheme which 
splits elongated cells. 



The so-called Orszag-Tang vortex ( Orszag Tang|[l979 ) is 
an often employed two-dimensional test for MHD codes. Be- 
sides being an excellent stability test it also examines how 
shocks interact with each other and produce a superson- 
ically turbulent system where the turbulence decays with 
time. The initial conditions as first described by |Picone fc| 
Dahlburg (1991) are defined by 



47r 
47r 



(20) 



V = (— sin(27r^), sin(27rx), 0) 
B = (-sin(27r2/),sin(47rx),0) 



in a box of extension < x, 2/ < 1- We consider test runs 
with a resolution of 600 x 600 moving mesh cells. Fig. |5] 



© 0000 RAS, MNRAS 000, 000-000 



Magnetohydrodynamics on an unstructured moving grid 7 

















^^^^^ 



Figure 5. Two-dimensional Orszag-Tang vortex test. Shown are linear density maps (pairs on the left) and magnetic energy density 
maps (pairs on the right) at times t = 0.25 (left) and t = 2.5 (right). For each pair, the left image shows the result of the AREPO 
code, while the right images displays the same simulation when the ATHENA code ( [Stone et al.|2Qd8| is used. The number of resolution 
elements is identical and set to 600 x 600 cells on both cases. The density ranges displayed in the color maps cover 0.06 to 0.5 (top left 
pair) and 0.1 to 0.4 (top right pair), respectively, and for the magnetic energy 0.0 to 0.4 (bottom left pair) and 0.0 to 0.3 (bottom right 
pair), respectively. 



compares the density and magnetic energy in our simula- 



tion with results obtained with the ATHENA code (Stone 
et al.|[2008 ), at two different times. In the earlier snapshot 
several shocks are just building up, while at the time of the 
second snapshot the simulation has already evolved for some 
time and become fully turbulent. A quantitative comparison 
at t = 0.5 is shown in Fig. |6] The agreement between the 
two different codes is in general very good and reassuring. 
Note that in contrast to the ATHENA code the magnetic ver- 
sion of AREPO does not keep the initial symmetry perfectly. 
This deviation from perfect symmetry is in part a result of 
the order in which the fluxes over different interfaces of a 
cell are applied. As this is more or less random in AREPO, 
small round-off asymmetries are introduced even for per- 
fectly symmetric conditions, and those tend to be amplified 
by means of the additional freedom allowed by the moving 
mesh. In addition, the ATHENA code shows slightly finer 
structures at later times. 



Due to the different underlying schemes, Voronoi grid in 
AREPO and Cartesian grid in ATHENA, there are also differ- 
ences in the computational resources the two codes require. 
In the Orszag-Tang vortex problem the AREPO consumes 
about three times the amount of memory the ATHENA code 
needs. This is mainly associated with the construction of the 
Voronoi grid which has to be done in each timestep. In ad- 
dition, ATHENA is up to six times faster than AREPO. This 
is caused by a combination of two effects. The timesteps 
in the AREPO code are smaller by a factor of about three 



due to stricter timestep criterions[^ and a single timestep 
takes twice as long for AREPO than for ATHENA on av- 
erage. The higher timestep cost is mainly caused by the 
construction of the Voronoi grid which takes about 40% of 
all CPU time, and the gradient and flux calculations which 
are slowed down considerably by cache misses in AREPO 
because the grid cells are stored unordered in memory. In 
addition, AREPO needs to solve more Riemann problems 
per cell. We note, however, that this difference is partic- 
ularly large in this application since the problem is very 
simple and dominated by magnetohydrodynamics only. For 
more complex applications that are often dominated by the 
evaluation of self-gravity this difference becomes much less 
important. 



4 DRIVEN SUPERSONIC MHD 
TURBULENCE 

Arguably one of the most challenging test problems for 
an MHD-code is driven supersonic MHD-turbulence, which 
we consider in this section. To this end, we create a pe- 
riodic box of unit size containing an isothermal gas with 
constant density. The initially uniform magnetic field has 
a strength of B = 0.1^3. In our dimensionless system of 
units, the volume, density and speed of sound are given by 
y = p = Cs = 1, respectively. To drive the turbulence, an 
external force field is applied to the fluid which we imple- 
ment as in [Price Federrath (2010). The force field is setup 



^ Using a less conservative timestep criterion can eliminate this 
contribution to the difference 
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Figure 6. Slice of the density at y = 0.33 and t = 0.5 of the 
Orszag-Tang vortex. Blue and red symbols show the results for 
the AREPO and the ATHENA code. 



Figure 7. The solid (dashed) lines show the sonic (Alfvenic) 
Mach number as a function of time. Both are shown for two sim- 
ulations with a resolution of 64^ cells (red) and 128"^ cells (black). 



in Fourier space and only contains power in a small range 
of low frequency modes between /cmin and /cmax- The am- 
plitude of the force field is given by a paraboloid with a 
maximum at kc = (/cmax — A^min)/2. The individual phases 
of the modes are drawn from Ornstein-Uhlenbeck processes. 
The stirring force field is then obtained grid by an inverse 
Fourier transformation at each cell center. 

The amplitude of the power spectrum of the driving 
forces is chosen to lead to an rms Mach number of about 
Ms ^5.5 once stationary turbulence has built up. The time 
evolution of this quantity, together wi th that of the approxi- 
mate Alfvenic Mach number Ma = \j2{pv?) / (B^) is shown 
in Figure |7| At the beginning, the Alfvenic Mach number 
increases rapidly, because the fluid velocity increases faster 
than the magnetic field. Later it decreases again as the mag- 
netic field catches up. We note that the sonic mach numbers 
of our two simulations with a resolution of 64^ and 128^ 
cells agree well with each other for an identical driving, as 
expected. However, the Alfvenic Mach number of the high 
resolution run is slightly smaller on average. The expected 
dynamical timescale in our setup is td — L/(2Ms), and we 
find that after about three dynamical time scales the tur- 
bulence is fully established. For the following analysis we 
therefore only consider outputs after t — Atd- 

To calculate power spectra of the turbulent velocity 
fields, we interpolate the Voronoi cells to a uniform Carte- 
sian grid using the Voronoi interpolation, or in other words, 
we sample the reconstructed velocity field defined on the 
Voronoi mesh at the coordinates of a fine Cartesian grid. To 
compare kinetic and magnetic power spectra, we also calcu- 
late the kinetic power spectra for a density- weighted velocity 
field, u = y^v. The magnetic power spectrum is obtained 
directly from the magnetic field. In this way, the integral 
over the power spectra is equivalent to the total kinetic and 
magnetic energy, respectively. 

The resulting magnetic and velocity power spectra are 
shown in Figure [S] They are averaged over 70 snapshots 
ranging from t — 4.4 td to 22.0 td- The velocity spectra at low 
wave numbers are largely dominated by the driving mech- 
anism. Towards smaller scales a power-law for the kinetic 
energy is found, followed by a dissipative regime when the 




Figure 8. Power spectra of driven MHD turbulence. Shown are 
the total magnetic power spectrum (solid lines), the divergence 
part of the magnetic power spectrum (dashed lines) and the ve- 
locity power spectrum (dotted lines). Spectra are shown for two 
runs with a resolution of 64^ cells (red) and 128^ cells (black), 
respectively. 



Nyquist frequency is approached. The inert ial range of the 
turbulence is however quite small at this comparatively low 
resolution. This probably also explains why no clear power 
law region for the magnetic power spectrum is seen. Note 
that in this setup the magnetic energy is about 1.5 orders 
of magnitude lower than the kinetic energy. The shape of 
the power spectra is similar to those found by |Kritsuk et al.| 
(2011) in a systematic comparison of different MHD codes, 
suggesting that our MHD implementation is consistent with 
the results found there. This is particularly encouraging as 
we here use a fully dynamic mesh, for which it is quite de- 
manding to control V • B errors under conditions of highly 
supersonic isothermal turbulence. 

We examine the latter point quantitatively by mea- 
suring the divergence part of the magnetic power spectra, 
obtained through the use of a Helmholtz decomposition in 
Fourier space. In the higher resolution run (128^), the V-B 
error is smaller by roughly an order of magnitude compared 
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with the low resolution run (64^), and in both cases, the 
error is very small compared to the actual strength of the 
magnetic field. 



5 COLLAPSE OF A MAGNETIZED CLOUD 

Finally, in this section we apply the code to the collapse of 
a magnetized sphere, which can be viewed as an idealized 
astrophysical problem that is of direct relevance to star for- 
mation. We choose initial conditions similar to the study of 
Hennebelle &: Fromang| |2008 ) . A rigidly rotating homoge- 
nous sphere with a radius of Ro = 0.015 pc and a total 
mass of one solar mass is embedded in an atmosphere that 
is 100 times less dense, with a small transition region at the 
boundary. With an initial density of 4.8 x 10~^^g cm~^ the 
freefall time is 3 x 10^ years. The sphere rotates rigidly with 
a period of 4.7 x 10^ years, equivalent to a ratio of rota- 
tional over gravitational energy equal to 0.045. The whole 
box in which the cloud is embedded is filled with a uniform 
magnetic field parallel to the rotation axis with a strength 
of 3.0 /xG, equivalent to a mass-to-flux over critical mass- 
to-flux ratio of 20. We adopt the same equation of state as 
Hennebelle k Fromang| ( [2008 ) given by 

P = px{c°f x[l + {p/p,f']°\ (21) 

"^^g cm~^. Our simula- 



with — 0.2 km s~^ and pc — 10~ 
tion box has a size of 0.06pc and inflow/outflow boundary 
conditions with an initial resolution of 128^ cells. We ap- 
ply a refinement criterion that splits a cell when its freefall 
timescale becomes smaller than ten times its sound crossing 
timescale. However, we limit the volume of a cell to be not 
smaller than 5 x 10~^^pc which is equivalent to an effective 
resolution of 16384^ cells. 

After 1.13tff the number of cells in the simulation in- 
creased from ^2.1 million in the beginning to ~ 2.2 mil- 
lion. In contrast to the initial setup, however, most of the 
cells are clustered in a very small region at the center of 
the cloud where most of the mass resides. The cloud at 
this time is shown in Fig. [9] As expected, a proto star has 
formed surrounded by an accretion disk with a radius of 
about 0.03 i?o- Magnetically powered outflows are launched 
along the z-axis with a velocity of 2. km s~^. The magnetic 
field parallel to the z-axis has been amplified by compres- 
sion beyond 10^ /xG close to the protostar. An azimuthal 
magnetic field has been generated during the collapse with 
a typical strength of 5/iG, up to 10^ /xG at the protostar, 
in good qualitative agreement with [Hennebelle &: Fromang] 
( 2008t . 

Fig. |9] also nicely shows how the adaptive mesh adapts 
to the internal structure of the fiuid in the simulation and 
is able to handle large density contrasts very well. 



interface, based on approximate solutions of Riemann prob- 
lems. To maintain the divergence constraint of the magnetic 
field we have implemented the divergence cleaning scheme 
proposed by Dedner et al. (2002|. In contrast to the con- 
strained transport approach, this method can be readily 
implemented for unstructured dynamic meshes that we use 
here. 

To our knowledge, our new MHD implementation is the 
first three-dimensional Lagrangian mesh code capable of fol- 
lowing magnetic fields. There already exist SPH implemen- 
tations of Lagrangian MHD (e.g. [Rosswog k Price] |2007 
|Dolag k Stasyszyn 2009 ), but they still suffer quite a bit 
from the inherent subsonic noise in SPH, necessitating rel- 
atively aggressive cleaning schemes, and from the relatively 
slow convergence rate of SPH (Springel |2010b| ). Our new 
scheme fairs much better in this respect, retaining the high 
accuracy for shocks and smooth fiows that is reached by 
ordinary fixed-mesh codes. In addition, our method dras- 
tically reduces advection errors in cases where large bulk 
velocities occur. Especially in this regime, our magnetic ver- 
sion of AREPO can be expected to surpass the accuracy of 
corresponding fixed mesh codes. 

We have verified the code's performance in a number 
of standard test problems, ranging from simple magnetic 
shock tubes to complicated interactions of multiple shock 
waves such as the Orszag-Tang vortex problem. In all our 
tests we found satisfactory agreement with published results 
from fixed-mesh MHD codes. We also applied our code to the 
collapse of a magnetized cloud of one solar mass under self- 
gravity, nicely reproducing results obtained by |Hennebelle| 



k Fromang (2008) with the RAMSES code. In addition, we 



have used our new code to simulate driven, highly supersonic 
isothermal magnetic turbulence in a periodic box. This has 
established that the new moving-mesh MHD code can han- 
dle this demanding problem robustly. The resulting power 
spectra are in good agreement to previous work considering 
that we have used a coarser resolution than in recent ded- 
icated MHD turbulence studies. In future work, it will be 
interesting to study the details of the statistics of the tur- 
bulence represented by the AREPO code, but this is beyond 
the scope of this paper and requires much larger production 
simulations. 

We conclude that our new MHD version of AREPO can 
be applied with confidence to a large range of astrophysical 
problems. Its Lagrangian properties should make it partic- 
ularly useful for studies of star formation, magnetic fields 
in galaxies and clusters of galaxies, and accretion disks. We 
also note that the code includes a powerful and accurate 
gravity solver, as well as modules for collisionless dynamics, 
making it also attractive for simulations of cosmic structure 
formation. 



6 CONCLUSIONS 

In this study, we have presented our implementation of 
ideal MHD in the moving-mesh code AREPO. The numer- 
ical scheme employs a fully adaptive Voronoi mesh which 
can freely move with the fiow. The dynamics is solved with 
a second-order accurate finite-volume scheme that employs 
a spatial reconstruction and a fiux calculation at each mesh 
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Figure 9. The collapsed cloud after 1.13 tff. Shown are the density (top row), inflow/outflow velocity (middle row), magnetic field in 
the z-direction (lower left plot) and in the azimuthal direction (lower right plot), respectively. 
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